#get samples
my_file="/Volumes/Data/copula_test/delay_exp14_1.txt";
X=scan(file=my_file);

#initialization
sampleSize=length(X);
outputSize=sampleSize*1000;
lagMax=10; #Max lag for autocor
Z<-seq(length=outputSize, 0, 0); #output standard gaussian
W<-seq(length=outputSize, 0, 0); #output

#add two extreme value
X[1]=-100000000000;
X[length(X)]=100000000000;

#map to standard gaussian
F=ecdf(X);
Fx=F(X);
Y<-seq(length=sampleSize-2, 0, 0);
for(i in seq(2, length(X)-1, 1))
{
	Y[i-1]=qnorm(Fx[i], mean=0, sd=1);
}
#remove first and last element of X
X=X[-length(X)];
X=X[-1];

#map X and Y
#sort, remove duplicated values in X and Y
Y_values=sort(unique(Y));
X_values=sort(unique(X));

#calculate autocor of the mapped standard gaussian
Ro=acf(Y, lag.max=lagMax, type="correlation", plot=FALSE);

#build covariance matrix (inversely), M means N-1
#build covariance matrix (inversely), M means N-1
Ro_N_N=seq(length=(lagMax+1)*(lagMax+1), 0, 0);
attr(Ro_N_N, "dim")=c(lagMax+1, lagMax+1);
for(i in seq(1,lagMax+1,1))
{
	for(j in seq(1, lagMax+1, 1))
	{
		Ro_N_N[i, j]=Ro$acf[abs(j-i)+1];
	}
}
Ro_1_2=Ro_N_N[c(1:lagMax),c((lagMax+1):(lagMax+1))];
attr(Ro_1_2, "dim")=c(lagMax,1);
Ro_2_1=Ro_N_N[c((lagMax+1):(lagMax+1)), c(1:lagMax)];
attr(Ro_2_1, "dim")=c(1,lagMax);
Ro_2_2=Ro_N_N[c((lagMax+1):(lagMax+1)), c((lagMax+1):(lagMax+1))];
Ro_1_1=Ro_N_N[c(1:lagMax), c(1:lagMax)];

#initialize output
mu_matrix=seq(length=(lagMax),0,0);
for(i in seq(1, lagMax, 1))
{
	Z[i]=Y[i];
	mu_matrix[i]=Y[i];
}
attr(mu_matrix, "dim")=c(lagMax,1);

temp1=Ro_2_1 %*% (solve(Ro_1_1));
new_var=Ro_2_2-(Ro_2_1 %*% (solve(Ro_1_1)) %*% Ro_1_2);
var_error=0;
if(new_var<=0)
{
	var_error=1;
}
new_ro=sqrt(new_var);

#generate Z from new mu and ro
for(i in seq(lagMax+1, outputSize, 1))
{
	new_mu=temp1 %*% mu_matrix;
	Z[i]=rnorm(1, mean=new_mu, sd=new_ro);
	
	for(j in seq(1, lagMax, 1))
	{
		mu_matrix[j,1]=Z[i-lagMax+j];
	}
}

#inverse back to W (actual output similar to X)
for(i in seq(1, outputSize, 1))
{
	found=0;
	for(j in seq(1, length(Y_values), 1))
	{
		if(Z[i]<=Y_values[j])
		{
			W[i]=X_values[j];
			found=1;
			break;
		}
	}
	
	if(found==0)
	{
		W[i]=X_values[length(X_values)];
	}
}

#plot cdf
plot(ecdf(X), col="red");
plot(ecdf(W), col="blue", add=TRUE);

quartz();
plot(ecdf(Y), col="red");
plot(ecdf(Z), col="blue", add=TRUE);

#calculate autocor
X_acf=acf(X, lag.max=lagMax, type="correlation", plot=FALSE);
W_acf=acf(W, lag.max=lagMax, type="correlation", plot=FALSE);
Y_acf=acf(Y, lag.max=lagMax, type="correlation", plot=FALSE);
Z_acf=acf(Z, lag.max=lagMax, type="correlation", plot=FALSE);

print(Y_acf);
print(Z_acf);
print(X_acf);
print(W_acf);